/* 
 * File:   surfalign_grdi3d.h
 * Author: boyan
 *
 * Created on April 9, 2009, 4:31 PM
 */

#ifndef _SURFALIGN_GRID3D_H
#define	_SURFALIGN_GRID3D_H

#include "AlignDomain.h"
#include "grid3d.h"
#include "LocalSurfaceCalculator.h"

#include <vector>

#include <MeshInterface.hpp>
#include <MsqError.hpp>

using Mesquite::Mesh;
using Mesquite::Vector3D;
using Mesquite::MsqError;

struct SurfFace;
struct SurfEdge;
struct SurfNode;
struct SurfNode1;

class SurfAlign_grid3d : public AlignDomain
{
public:
    explicit SurfAlign_grid3d(grid3d_t *g);
    SurfAlign_grid3d(grid3d_t *g, std::vector<LevelsetCode> &ls);

    virtual ~SurfAlign_grid3d()
    {
    }

    // Inherited from AlignDomain
    virtual LevelsetCode levelset_code(Mesquite::Mesh::EntityHandle entity_handle,
            const Mesquite::Vector3D &coordinate) const;

    // Inherited from MeshDomain
    virtual void snap_to(Mesh::EntityHandle entity_handle,
            Vector3D &coordinate) const;
    virtual void normal_at(Mesh::EntityHandle entity_handle,
            Vector3D &coordinate) const;
    virtual void normal_at(const Mesh::EntityHandle* handles,
            Vector3D coordinates[], unsigned count, MsqError& err) const;
    virtual void closest_point(Mesh::EntityHandle handle,
            const Vector3D& position, Vector3D& closest,
            Vector3D& normal, MsqError& err) const;

    virtual void domain_DoF(const Mesh::EntityHandle* handle_array,
            unsigned short* dof_array, size_t num_handles, MsqError& err) const
    {
        for(size_t i = 0; i < num_handles; i++)
            dof_array[i] = 2;
    }


    virtual void write_vtk(const char* filename) const;
    virtual void write_vtk(Mesh *mesh, const char* filename,
            bool write_codes = false);

    virtual void get_codes(const Mesh::EntityHandle *handles,
            unsigned codes[], size_t count)
    {
        for(size_t i = 0; i < count; i++) codes[i] = 0;
    }

protected:
    struct GridElem;
    grid3d_t *grid;
    std::vector<GridElem> gridElems;


    std::vector<int> inside_elems;
    std::vector<int> nodes_inv;
    std::vector<int> edges_inv;
    std::vector<int> faces_inv;
    std::vector<SurfFace> surf_faces;
    std::vector<SurfEdge> surf_edges;
    std::vector<SurfNode> surf_nodes;

    // y(x) = surf_nodes[x].id is index in grid->nodes
    // x(y) = nodes_inv[y] is index in surf_nodes
    // nodes_inv[surf_nodes[x].id] = x,  for 0 <= x < surf_nodes.size()
    // surf_nodes[nodes_inv[y]].id = y,  for 0 <= y < grid->NN if y is node in surf_nodes
    //                             = undefined, if y is not a node in surf_nodes
    //     if y is not in surf_nodes, then nodes_inv[y] = -1, which is garbage
    // Similar situation with egdes and faces.

    void initialize();
    void build_faces();
    void build_nodes();
    void build_edges();
    void __compute_ball(SurfNode &X,
            const std::vector<size_t> *edges = NULL,
            const std::vector<size_t> *faces = NULL) const;

    void __compute_ball1(SurfNode1 &X,
            const std::vector<size_t> *edges = NULL,
            const std::vector<size_t> *faces = NULL) const;

    friend bool __radius_centre_of_ball(SurfNode1 &X, std::vector<Vector3D> &pts);

    virtual size_t find_closest_node(Mesh::EntityHandle handle,
            const Vector3D position) const;
    virtual bool isPointInside(const Vector3D &pos) const;
    // ************************************************************************
    std::vector< std::list<int> > AdjacentNodes;
    std::vector< std::list<int> > AdjacentNodes2;
    std::vector< std::list<int> > AdjacentNodes3;

    void createAdjacency();
};

struct SurfFace
{
    int id;
    bool swap_normal;
    static grid3d_t *grid;

    explicit SurfFace(int fid)
    : id(fid), swap_normal(false)
    {
    };

    Vector3D getFaceNormal(double *Area = NULL) const
    {
        const int *my_nodes = nodes();
        Vector3D n[3];
        for(int i = 0; i < 3; i++)
            n[i].set(grid->nodes[my_nodes[i]].coords);
        Vector3D normal = // cross product; the order matters
                (n[1] - n[0]) * (n[2] - n[0]);
        double A = normal.length();
        if(Area) *Area = A;
        normal /= A; // normal.normalize();
        if(swap_normal)
            return -normal;
        else
            return normal;
    }

    bool has_node(int nid) const
    {
        const int *my_nodes = nodes();
        return (my_nodes[0] == nid) || (my_nodes[1] == nid) || (my_nodes[2] == nid);
    };

    int surf_nodes(int i, const std::vector<int> &inv) const
    {
        assert((int) inv.size() == grid->NN);
        return inv[nodes(i)];
    }

    int surf_edges(int i, const std::vector<int> &inv) const
    {
        assert((int) inv.size() == grid->NE);
        return inv[edges(i)];
    }

    const int *nodes() const
    {
        return grid->faces[id].nodes;
    };

    const int *edges() const
    {
        return grid->faces[id].edges;
    };

    const int *elems() const
    {
        return grid->faces[id].elems;
    };

    int nodes(int i) const
    {
        return grid->faces[id].nodes[i];
    };

    int edges(int i) const
    {
        return grid->faces[id].edges[i];
    };

    int elems(int i) const
    {
        return grid->faces[id].elems[i];
    };

};

struct SurfEdge
{
    int id;
    static grid3d_t *grid;

    explicit SurfEdge(int eid)
    : id(eid)
    {
    };

    bool has_node(int nid) const
    {
        const int *my_nodes = nodes();
        return (my_nodes[0] == nid) || (my_nodes[1] == nid);
    };

    bool has_surf_node(int nid, const std::vector<int> &inv) const
    {
        const int *my_nodes = nodes();
        return (inv[my_nodes[0]] == nid) || (inv[my_nodes[1]] == nid);
    }

    int other_node(int nid) const
    {
        const int *my_nodes = nodes();
        return my_nodes[0] + my_nodes[1] - nid;
    }

    int other_surf_node(int nid, const std::vector<int> &inv) const
    {
        const int *my_nodes = nodes();
        return inv[my_nodes[0]] + inv[my_nodes[1]] - nid;
    }

    const int *nodes() const
    {
        return grid->edges[id].nodes;
    };

    int nodes(int i) const
    {
        return grid->edges[id].nodes[i];
    };

    int surf_nodes(int i, const std::vector<int> &inv) const
    {
        assert((int) inv.size() == grid->NN);
        return inv[nodes(i)];
    }

};

struct SurfNode
{
    int id;
    static grid3d_t *grid;
    PointSetSurfaceCalculator * calc;

    explicit SurfNode(int nid) : id(nid), calc(0)
    {
    }

    SurfNode(const SurfNode & to_copy) : id(to_copy.id), calc(NULL)
    {
        // do not merely copy the pointer, reather create a new copy of the instance
        if(to_copy.calc)
            setCalc(*(to_copy.calc));
    }

    SurfNode &operator=(const SurfNode &to_copy)
    {
        id = to_copy.id;
        if(to_copy.calc)
            setCalc(*to_copy.calc);
        else
            calc = NULL;
        return *this;
    }

    virtual ~SurfNode()
    {
        if(calc)
            delete calc;
    }

    void setCalc(const PointSetSurfaceCalculator & like_me)
    {
        if(calc)
            delete calc;
        calc = like_me.NewCopy();
        calc->setFiducialPoint(toVector3D());
    };

    void flipNormal()
    {
        if(calc)
            calc->flipNormal();
    }

    const double *toArray() const
    {
        return grid->nodes[id].coords;
    }

    Vector3D toVector3D() const
    {
        return Vector3D(grid->nodes[id].coords);
    }

    Vector3D GetSurfNormal() const
    {
        return calc->getNormal();
    }

    Vector3D GetSurfNormal(const Vector3D & pt) const
    {
        return calc->getNormal(pt);
    }

    Vector3D GetProjection(const Vector3D & pt) const
    {
        return calc->getProjection(pt);
    }
};

struct SurfNode1
{
public:
    int id;
    double Kappa; // curvature == 1/Radius;
    Vector3D BallCenter;
    // size_t EdgesNum;
    signed int NormalSign;
    static grid3d_t *grid;

    explicit SurfNode1(int nid)
    : id(nid),
    Kappa(0.0), BallCenter(0.0, 0.0, 0.0), // EdgesNum(0),
    NormalSign(0)
    {
    }

    SurfNode1(const SurfNode1 & rhs)
    : id(rhs.id),
    Kappa(rhs.Kappa), BallCenter(rhs.BallCenter), // EdgesNum(rhs.EdgesNum),
    NormalSign(rhs.NormalSign)
    {
    }

    const double *toArray() const
    {
        return grid->nodes[id].coords;
    }

    Vector3D toVector3D() const
    {
        return Vector3D(grid->nodes[id].coords);
    }

    SurfNode1 & operator=(const SurfNode1 & rhs)
    {
        id = rhs.id;
        Kappa = rhs.Kappa;
        BallCenter = rhs.BallCenter;
        // EdgesNum = rhs.EdgesNum;
        NormalSign = rhs.NormalSign;
        return *this;
    }

    Vector3D GetSurfNormal() const
    {
        return (NormalSign < 0) ? -GetBallNormal() : GetBallNormal();
    }

    Vector3D GetSurfNormal(const Vector3D & pt) const
    {
        if(Kappa > 0)
            return (NormalSign < 0) ? -GetBallNormal(pt) : GetBallNormal(pt);
        else
            return GetSurfNormal();
    }

    Vector3D GetBallNormal() const
    {
        Vector3D N = toVector3D() - BallCenter;
        N.normalize();
        assert(std::isfinite(N.x()));
        return N;
    }

    Vector3D GetBallNormal(const Vector3D & pt) const
    {
        Vector3D N = pt - BallCenter;
        N.normalize();
        if(!std::isfinite(N.x()))
            return GetBallNormal();
        else
            return N;
    }

    Vector3D GetProjection(const Vector3D & pt) const
    {
        if(Kappa > 0)
        {
            // let's see if point is above or below equator if *this is the North Pole
            if((pt - BallCenter) % (toVector3D() - BallCenter) > 0)
                return BallCenter + GetBallNormal(pt) / Kappa;
            else
                return BallCenter - GetBallNormal(pt) / Kappa;
        }
        else
        {
            Vector3D N = GetBallNormal();
            double alpha = (pt - toVector3D()) % N;
            return pt - alpha*N;

            // plane projection reminder: Gram-Schmidt
            // plane with point N and normal n.
            // projection of point P => C
            // the formula is: C = P - alpha n
            // where alpha = <NP, n> / <n,n>, where vector NP = P - N
        }
    }

    bool IsNormalSignSet() const
    {
        return NormalSign != 0;
    };

};

struct SurfAlign_grid3d::GridElem
{
    int id;
    double box[2][3];
    static grid3d_t *grid;

    explicit GridElem(int Tid)
    : id(Tid)
    {
        elem_t *T = grid->elems + id;
        double *nds[4];
        for(int i = 0; i < 4; i++)
            nds[i] = grid->nodes[T->nodes[i]].coords;
        for(int i = 0; i < 3; i++)
        {
            box[0][i] = fmin(fmin(nds[0][i], nds[1][i]), fmin(nds[2][i], nds[3][i]));
            box[1][i] = fmax(fmax(nds[0][i], nds[1][i]), fmax(nds[2][i], nds[3][i]));
        }
    };

    GridElem(const SurfAlign_grid3d::GridElem & to_copy)
    {
        id = to_copy.id;
        for(int i = 0; i < 2; i++)
            for(int j = 0; j < 3; j++)
                box[i][j] = to_copy.box[i][j];
    }

    inline bool contains_point(const Vector3D & pt) const
    {
        // is pt it in box?
        const double *p = pt.to_array();
        for(int i = 0; i < 3; i++)
            if(p[i] < box[0][i] || p[i] > box[1][i])
                return false;
        // pt is in box.  is in it element?
        double t[4];
        barycentric_coords(pt, t);
        for(int i = 0; i < 4; i++)
            if(t[i] < -1e-6)
                return false;
        return true;
    }

    inline void barycentric_coords(const Vector3D &pt, double t[4]) const
    {
        int i, j, *fnds;
        double C, A, N[3];

        elem_t *T = grid->elems + id;

        t[3] = t[2] = t[1] = t[0] = 0.0;
        C = 1.0 / 3.0 / ComputeVolume(grid->nodes, T->nodes);
        for(i = 0; i < 4; i++)
        {
            fnds = elem_face_nodes(T->nodes, i);
            A = ComputeArea(grid->nodes, fnds, N);
            for(j = 0; j < 3; j++)
                t[i] += (grid->nodes[fnds[0]].coords[j] - pt[j]) * N[j];
            t[i] *= A * C;
        }
    }
};


#endif	/* _SURFALIGN_GRID3D_H */

